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1  Research  Objective 

Understanding  the  failure  of  metals  at  high  strain  rate  is  of  utmost  importance  in  the  design  of 
a  broad  range  of  Army  systems.  Numerical  methods  offer  the  ability  to  analyse  such  complex 
physics  and  aid  the  design  of  structural  systems.  The  objective  of  this  research  was  to  develop 
reliable  finite  element  models  for  high  strain  rate  failure  modeling,  simultaneously  incorporating 
shear  banding  and  fracture.  To  our  knowledge,  there  are  no  such  models  that  combine  both  shear 
bands  and  fracture  under  one  unified  framework,  which  represents  a  significant  technical  barrier  to 
detailed  failure  prediction  and  design  of  material  systems  under  high  rate  loading. 

To  this  end,  we  have  developed  a  computational  framework  that  includes  monolithic,  mixed 
finite  element  formulation  which  accounts  for  both  phenomena  in  a  regularized  continuum  sense 
leading  to  mesh-insensitive  results. 

Shear  bands  are  driven  by  shear  heating  caused  by  isochoric  inelastic  deformations  and  are 
regularized  with  thermal  diffusion.  Cracks  propagate  when  the  internal  strain  energy  reaches  a 
critical  value  and  are  regularized  by  using  a  phase  field  model,  which  models  a  crack  as  a  diffusive 
entity.  Both  shear  bands  and  cracks  have  been  identified  as  failure  modes  in  metals  at  high  strain 
rate.  We  therefore  contend  that  such  a  model  will  vastly  improve  quantitative  failure  prediction 
for  structural  components  loaded  at  high  strain  rates  and  will  support  Army  research  with  novel 
unified  framework  for  shear  bands  and  cracks. 

Several  Army  applications  in  lethality  and  protection  can  be  enhanced  by  the  proposed  model. 
Quantitative  prediction  of  the  energy  absorption  capability  of  protection  systems,  vehicle  struc¬ 
tures,  and  welded  joints  is  extremely  important  for  evaluation  of  the  survivability  under  extreme 
loading  conditions.  Kinetic  energy  penetrators  are  another  area  of  interest,  where  enhanced  lethal¬ 
ity  can  be  achieved  through  quantitative  prediction  of  self  sharpening  behavior.  It  has  been  ob¬ 
served  that  both  fracture  and  shear  bands  are  the  main  failure  mechanisms  under  high  rate  loading. 
Thus,  the  use  of  numerical  models  which  exclude  either  failure  mechanism  will  tend  to  over  es¬ 
timate  the  energy  absorption  capacity.  The  proposed  model  seeks  to  rectify  this  deficiency.  A 
particular  ARL  application  of  interest  is  the  failure  analysis  of  friction  stir  welded  (FSW)  alu¬ 
minum  joints  for  use  in  Army  land  vehicles.  This  further  introduces  material  inhomogeneity  and 
phase  interface  issues  between  various  weld  zone  microstructures,  all  of  which  requires  detailed 
dedicated  analysis. 

2  Introduction:  Dynamic  fracture  and  Shear  Bands 

2.1  Literature  review 

Metals  and  alloys  subjected  to  dynamic  loads  such  as  impact  or  blast,  may  fail  due  to  brittle  and/or 
ductile  fracture,  shear  banding,  or  combination  of  the  two  depending  on  factors  such  as  material 
properties,  loading  rate  and  specimen  geometry. 

Fracture  and  shear  banding  are  two  different  phenomenon  in  solid  mechanics  with  distinct 
spatial  and  temporal  scales.  Shear  bands  are  zones  of  highly  localized  plastic  deformation,  which 
in  metals  are  usually  observed  in  high  strain  rate  loadings,  while  brittle/ductile  cracks  are  material 
discontinuities  due  to  cleavage  and/or  void  coalescence. 

Shear  bands  have  also  been  observed  in  materials  with  non-associative  flow  laws  like  porous  or 
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granular  materials,  polymers,  rocks  and  soils  [1,2].  High  strain  rate  loading  amplifies  strain  rate 
hardening  effects  observed  in  many  structural  materials,  which  has  a  stabilizing  effect.  In  contrast, 
thermal  softening  destabilizes  the  material  as  heat  produced  by  inelastic  deformation  lowers  the 
material’s  flow  stress.  The  reduced  flow  stress  then  leads  to  more  intense  plastic  deformation  where 
an  instability  may  occur.  The  influence  of  strain  hardening,  strain  rate  hardening  and  thermal  soft¬ 
ening  produces  a  complex  equilibrium  that  depends  strongly  on  the  loading  rate  and  the  material 
properties.  Hence,  thermal  softening  is  well  accepted  as  the  main  mechanism  that  leads  to  shear 
band  formation[3-8].  However,  recent  studies  suggest  that  dynamic  recrystallization  (DRX)  may 
be  another  potential  mechanism  that  can  trigger  shear  band  failure  [9-1 1],  but  it  is  not  considered 
in  the  current  work. 

While  shear  banding  is  often  the  main  driver  of  failure  by  ductile  fracture[12],  it  is  also  consid¬ 
ered  a  failure  mode  in  its  own  right,  since  the  thermal  softening  leads  to  profound  and  rapid  loss 
of  load  carrying  capability  [3,  13,  14]. 

Experimentally  derived  material  models  for  these  loading  regimes  describe  plastic  flow  as  be¬ 
ing  dependent  on  temperature,  strain  rate,  and  a  hardening  parameter  [7].  While  several  models  are 
available  for  modeling  shear  bands,  most  are  similar  in  that  increasing  temperature  (due  to  plastic 
work)  has  a  softening  effect  which  cause  plastic  flow  to  occur  more  readily,  while  increases  in 
strain  rate  and  the  hardening  parameter  have  a  hardening  effect.  Following  the  experimental  work 
of  Marchand  and  Duffy  [15],  shear  bands  develop  in  three  stages  as  shown  in  Figure  1.  In  Stage 
I,  before  localization,  a  homogeneous  distribution  of  plastic  strain  is  observed.  Stage  II  begins 
when  the  thermal  softening  effect  dominates  the  strain  and  strain  rate  hardening  effects,  resulting 
in  strain  softening,  and  thus  mild  strain  localization.  Stage  III  is  marked  by  severe  localization  and 
rapid  softening,  a  phenomena  termed  stress  collapse,  which  indicates  a  sudden  and  large  drop  in 
the  material’s  load  bearing  capability  [5]. 


Figure  1 :  Stages  of  Shear  band  deformation.  I-homogeneous  solution,  Il-onset  of  localization. 
Ill-stress  collapse. [15] 


Since  shear  bands  are  typically  driven  by  shear  heating,  a  Taylor-Quinney(TQ)  parameter  is 
typically  considered.  This  coefficient  determines  the  fraction  of  inelastic  work  that  is  converted 
to  heat  and  in  most  studies  in  the  literature  has  been  assumed  to  be  constant  with  a  value  of  0.9. 
Nonetheless,  recent  studies  indicate  that  the  TQ  coefficient  is  not  constant  and  may  depend  on 
internal  variables  such  as  strain  rate,  plastic  strain  and  temperature,  among  others.  These  stud¬ 
ies  assume  that  some  portion  of  the  plastic  work  goes  into  cold  work,  which  is  energy  stored  in 
defects  within  the  metal’s  lattice  [16,  17].  Following  these  studies  we  will  consider  in  this  work 
a  generic  non-constant  strain-rate  dependent  Taylor  Quinney  coefficient.  Note  that  the  functional 
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non-linearity  of  the  TQ  coefficient  is  an  active  field  of  research  and  proper  satisfactory  consensus 
is  yet  to  be  found  [16-19]. 

The  interplay  between  inelastic  deformation  associated  with  shear  bands  and  subsequent  duc¬ 
tile  fracture  due  to  coalescence  and  growth  of  voids  has  been  another  topic  of  intense  research. 
Brittle-ductile  failure  transition  has  been  observed  in  the  impact  of  notched  plates  by  [20-22].  It 
was  found  that  failure  by  shear  banding  occurred  above  a  critical  impact  velocity,  while  brittle  frac¬ 
ture  resulted  from  lower  impact  velocities.  While  this  counterintuitive  phenomenon  is  triggered  by 
strain  rate  effects,  the  more  well  known  (reverse)  ductile-brittle  failure  transition  occurs  due  to 
decrease  in  ambient  temperature [23].  Furthermore,  it  is  interesting  to  note  that  at  a  transition  ve¬ 
locity,  both  failure  modes  could  be  present  at  the  same  time  and  any  reliable  model  should  attempt 
to  capture  both  modes. 

Another  important  example  is  the  impact  and  penetration  of  a  metal  plate  by  a  projectile.  In 
this  experiment,  shear  bands  have  been  identified  propagating  parallel  to  the  impact  direction,  and 
cracks  radiate  perpendicular  to  the  impact  direction  [24].  Again,  both  failure  modes  are  present 
at  the  same  time.  These  failure  modes  present  in  dynamic  fracture  problems  involve  different 
phenomenon  with  distinct  spatial  and  temporal  scales.  Hence  to  model  and  characterize  metals 
failure  under  dynamic  loading  conditions,  it  is  thus  crucial  to  account  for  all  failure  modes,  since 
exclusion  of  any  one  of  them  neglects  important  underlying  physics. 

To  this  end,  recent  significant  efforts  have  been  made  to  combine  the  effects  of  plasticity  and 
ductile  fracture  through  the  phase-field  method  [25,  26].  Such  a  model  has  been  recently  proposed 
by  McAuliffe  and  Waisman  [27,  28],  where  a  unified  single  set  of  PDFs  is  used  to  capture  brittle 
and/or  ductile  fracture  and  shear  bands.  In  that  approach,  fracture  is  modeled  with  a  novel  phase 
field  formulation  and  is  coupled  to  the  set  of  equations  that  describe  shear  bands  [3,  13,  29-31],  as 
illustrated  schematically  in  Figure  2. 


(a)  Experimental  Results  [24] 


(b)  Phase-field  with  shear  band  modeling[27] 


Figure  2:  Projectile  penetrating  a  metal  plate  showing  both  fracture  and  shear  localization.  Fig¬ 
ure  2a  shows  the  experimental  results  of  Nickodemus  et  al.  [24]  where  the  plate  is  failing  due  to 
shear  deformation  in  the  direction  parallel  to  the  penetration  and  cracks  nucleate  and  propagate 
perpendicularly  to  the  projectile  direction.  In  Figure  2b,  an  illustration  of  the  unified  model  pre¬ 
sented  in  McAuliffe  and  Waisman  [27]  is  depicted,  where  the  phase  field  approximation  to  a  crack 
is  shown. 


The  phase-field  technique  is  a  general  methodology  that  has  been  extensively  used  to  model 
phase  transformations  in  materials,  particularly  when  dealing  with  problems  with  moving  bound¬ 
aries.  These  phases  can  be  associated  with  several  distinct  physical  phenomena,  that  include  so¬ 
lidification  [32-34],  solid-state  phase  transformation  [35,  36],  martensitic  transformation  [37,  38], 
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dislocation  dynamics  [39],  grain  growth  [40, 41],  among  others.  In  this  work  we  use  the  phase-field 
formulation  within  the  scope  of  fracture  mechanics,  which  has  been  extensively  studied  [42-54]. 

As  depicted  in  figure  3,  in  the  phase-field  formulation  a  crack  is  no  longer  modeled  as  a 
strong  discontinuity  in  the  domain  but  is  instead  approximated  by  a  continuous  surface  density 
function[25,  26,  42-46,  48,  49].  The  extent  of  the  damage  to  the  material  is  characterized  by  the 
phase  field  parameter  c,  which  assumes  the  value  of  1  in  the  fully  fractured  phase  and  0  in  the 
undamaged  phase,  and  the  width  of  the  diffusive  region  is  determined  by  the  length  scale  param¬ 
eter  Iq  [44,  55].  This  formulation  has  the  so-called  property  of  F-convergence  in  Iq,  which  means 
that  the  results  will  converge  to  standard  discrete  fracture  mechanics  results  when  Iq  — ?•  0[50,  56]. 
The  value  of  /q  is  typically  set  empirically  and  smaller  values  require  a  larger  density  of  elements 
near  the  crack.  However,  some  authors  have  argued  that  /q  should  be  treated  instead  as  a  material 
parameter  and  calibrated  based  on  experimental  results  [57]. 


Figure  3:  Schematic  depiction  of  a  solid  with  a  discontinuity  F  (crack).  In  the  phase-field 
formulation  the  crack  is  simulated  by  the  field  c  where  a  black  color  corresponds  to  a  fully  damaged 
material  (c  =  1)  and  a  white  color  corresponds  to  a  fully  intact  material  (c  =  0).  The  Iq  parameter 
controls  the  width  of  the  process  zone. 

One  of  the  most  important  features  that  control  the  behavior  of  the  phase  field  method,  is  the 
so-called  degradation  function.  This  function  is  used  to  degrade  a  chosen  component  of  the  elastic 
strain  energy  (e.g.  the  energy  associated  with  tension)  and  is  typically  associated  with  the  current 
state  of  damage  (i.e.  the  phase-field  parameter). 

The  degradation  function,  denoted  by  m(c),  can  be  chosen  arbitrarily  as  long  as  it  satisfies  the 
conditions  necessary  to  observe  Gamma  convergence  [56]  and  boundedness  of  the  fracture  force 
[48].  In  the  bulk  of  the  work  done  on  phase  field  methods,  a  quadratic  degradation  function  is  typ¬ 
ically  used  as  it  is  the  only  one  for  which  F-convergence  has  been  demonstrated  [50].  However,  in 
recent  years  other  degradation  functions  have  been  developed  [57].  For  example  a  cubic  degrada¬ 
tion  function  has  been  suggested  by  Borden  [58]  and  was  shown  to  provide  some  advantages  like 
closer  resemblance  to  an  elastic-brittle  behavior  as  expected  and  obtained  by  traditional  methods 
and  reduction  of  the  mild  degradation  that  occurs  away  from  the  crack. 

2.2  Modeling  aspects 

Given  the  short  time  scales  associated  with  the  formation  of  shear  bands,  adiabatic  conditions  have 
been  commonly  used  in  order  to  simplify  the  modeling  of  the  problem.  However,  this  assumption 
also  leads  to  mesh-dependent  results  [30,  59-61].  By  considering  the  effect  of  thermal  diffusion. 
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an  intrinsic  length  scale  is  introduced  that  functions  as  a  localization  limiter,  and  leads  to  mesh 
insensitive  results  [4,  62,  63].  Other  regularization  techniques  reported  in  the  literature  include 
strain  gradient  theories  [64-66],  which  have  been  used  in  the  context  of  shear  bands  in  [67-69] 
and  nonlocal  methods  [70-72].  Mesh  alignment  is  another  form  of  sensitivity  that  can  be  improved 
using  mesh-free  formulations [73,  74]  or  Isogeometric  analysis  [75]. 

The  phase-field  approach  results  in  a  mathematical  formulation  for  a  diffused  crack  that  is 
closely  related  to  gradient  damage  mechanics  [56,  76,  77],  although  it  has  been  argued  that  this 
correspondence  is  somewhat  coincidental  since  the  approaches  focus  on  modeling  different  phe¬ 
nomena  [76,  78].  The  significant  difference  between  a  phase  field  approach  and  a  discrete  crack 
approach  is  in  the  representation  of  a  crack,  continuous  versus  discrete,  and  would  mostly  depend 
on  the  application  problem.  Nonetheless,  continuous  representation  of  cracks  has  the  advantage 
that  it  allows  the  functional  minimization  problem  to  be  solved  by  numerical  methods,  such  as  the 
Finite  Element  Method,  without  the  need  for  any  special  set  of  shape-functions  or  enrichments. 
Hence,  complex  crack  patterns  such  as  branching  and  coalescence  can  easily  be  captured  with  this 
approach.  Even  though  this  technique  has  been  initially  derived  to  tackle  brittle  fracture,  it  has 
also  been  employed  to  model  brittle,  quasi-brittle  and  ductile  fracture  combined  with  plasticity 
[27,  45,  79,  80]. 

In  summary,  two  length  scales  are  present  in  this  model:  a  strong  length  scale  Iq  dominated 
by  fracture  and  a  weak  length  scale  dominated  by  shear  bands.  These  length  scales  regularize  the 
PDE  system  and  have  been  shown  to  provide  reliable,  mesh  insensitive  results  [81,  82]. 


2.3  Stability  Analysis 

In  this  work  we  also  formulate  stability  conditions  for  a  dynamic  fracture  and  shear  band  local¬ 
ization  model[27,  28].  Most  of  the  analysis  is  done  using  the  linear  perturbation  method  but  other 
techniques  are  also  proposed  and  studied. 

To  motivate  the  stability  analysis,  consider  that  in  a  homogeneous  problem  (e.g.  a  bar  in  tension 
without  imperfections)  the  phase-field  formulation  has  two  distinct  behaviors.  Initially,  the  entire 
problem  behaves  homogeneously  with  initiation  of  damage  due  to  the  driving  term  of  the  phase- 
field  equation  [83].  After  a  certain  critical  point,  the  system  loses  stability  and  a  non-homogeneous 
deformation  is  possible,  eventually  leading  to  fracture.  Hence,  loss  of  stability  provides  important 
information  and  is  typically  associated  with  the  onset  of  excessive  deformation  or  damage  and 
decreased  reliability  of  the  results  in  numerical  simulations 

In  quasi-static  simulations  stability  can  be  affected  by  phenomena  like  bifurcation  due  to  non- 
associated  flow  law  [1,  84,  85]  and  will  limit  the  robustness  of  the  numerical  approach,  often 
requiring  more  advanced  stepping  procedures  (e.g.  the  arc-length  method  [86])  and  regularization 
techniques  to  prevent  mesh  dependency  in  localization  phenomena  [4,  61,  67,  87].  In  dynamic 
problems  stability  is  often  associated  with  the  appearance  of  a  non-homogeneous  type  of  solution, 
for  example  in  thermo-mechanical  shearband  localization  [3,  88,  89]. 

Some  examples  of  the  importance  of  stability  detection  can  be  found  in  the  work  of  Rabczuk 
et  al.  [90]  that  set  strong  displacement  discontinuities  with  cohesive  surfaces  where  local  instability 
is  detected.  Song  et  al.  [91]  introduce  additional  tangential  degrees  of  freedom  at  the  discontinuity 
using  phantom  nodes.  Belytschko  et  al.  [92]  inject  coarse  scale  discontinuities  based  on  isolating 
the  unstable  regions  at  the  fine  scale  into  a  localization  domain.  Tabarraei  et  al.  [93]  use  material 
instability  based  on  the  perturbation  method  to  inject  a  phantom  node  discontinuity  on  a  macro 
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level  model  based  on  the  cohesive  laws  obtained  from  a  microstructural  model.  Rabczuk  and 
Samaniego  [94]  use  instability  in  an  adaptive  mesh-free  method  for  discontinuous  modeling  of 
shear  bands.  Finally,  Belytschko  et  al.  [95]  set  strong  displacement  discontinuities  with  a  cohesive 
law  where  local  instability  is  detected  for  fracture. 

A  popular  technique  for  analyzing  the  stability  of  a  system  is  the  linear  perturbation  method 
which  has  been  applied  successfully  to  shear  band  problems[3,  13,  14,  63,  96-102]. 

We  note  that  another  popular  stability  approach  that  is  widely  used  in  the  literature  is  based 
on  an  eigenvalue  analysis  of  the  acoustic  tensor  [103-106].  While  the  null-space  criterion  of  the 
acoustic  tensor  has  been  used  to  study  the  stability  of  failure  processes  in  rate  independent  mate¬ 
rials,  it  has  been  suggested  in  the  literature  that  this  approach  is  not  well  suited  for  rate  dependent 
materials  (as  in  our  case)  as  the  equations  remain  elliptic  [107,  108].  Therefore  we  choose  not  to 
pursue  this  route  and  focus  on  the  linear  perturbation  method. 

However,  the  success  of  the  linear  perturbation  method  comes  at  a  significant  analytical  cost. 
New  cases  must  undergo  the  laborious  process  of  deriving  a  stability  condition  based  on  the  gov¬ 
erning  equations  and  often  the  extension  of  these  criteria  to  multi-dimensions  is  not  trivial.  Hence, 
in  this  work  we  also  propose  a  method  to  determine  the  onset  of  shear  band  instability  by  a  spe¬ 
cially  formulated  local  generalized  eigenvalue  analysis.  To  this  end,  the  semi-discrete  Jacobian  is 
decomposed  from  its  mixed  finite  element  formulation  following  the  derivation  in  [81].  Although 
similar  techniques  that  compute  the  eigenvalues  of  the  Jacobian  matrix  originated  from  the  finite 
element  discretization  have  been  employed  in  other  fields  (e.g.  stability  of  fluid  dynamics  [109] 
or  thermal  capillarity  instability  [110]),  this  methodology  has  not  been  employed  to  study  local 
instabilities  of  shear  bands  before.  It  is  shown  that  the  local  stability  criteria  based  on  the  pertur¬ 
bation  method  mentioned  above,  under  the  same  assumptions,  is  equivalent  to  the  local  eigenvalue 
analysis  [88]. 

Like  the  two  techniques  used  in  this  work  and  described  above,  the  stability  of  a  system  is 
traditionally  studied  from  the  perspective  of  the  modal  spectrum  of  that  system.  The  underly¬ 
ing  reasoning  is  that  whichever  mode  has  the  largest  growth  rate  corresponding  to  an  exponential 
growth  will  dominate  the  solution.  However,  for  non-normal  systems  (i.e.  systems  whose  eigen¬ 
vectors  are  non-orthogonal)  the  short  time  growth-rate  of  a  solution  might  not  correspond  to  the 
growth  rate  of  the  most  unstable  mode.  On  the  other  hand,  non-autonomous  systems  do  not  have 
a  constant  spectrum  because  of  the  time  dependency  of  their  operators.  Since  the  problem  being 
studied  is  neither  normal  nor  autonomous,  a  generalized  stability  analysis  is  of  interest  [111,  112]. 

There  are  several  different  concepts  that  are  important  for  a  generalized  stability  analysis  as 
shown  in  Farrell  and  loannou  [111,  112].  The  study  of  the  growth  rate  of  a  perturbation  tangent  to 
the  solution  is  of  particular  interest  since  it  is  computationally  cheap  and  can  easily  be  implemented 
on  complex  multiphysics  systems.  Additionally,  an  analysis  based  on  the  concept  of  the  Rayleigh 
Quotient  of  an  eigenvector  is  employed  and  is  shown  to  also  produce  a  good  approximation  of  the 
instability  point. 

Albeit  considerably  more  expensive,  the  local  spectral  analysis  posses  the  same  implementation 
simplicity  that  the  instantaneous  growth  rate  and  the  Rayleigh  Quotient  have,  i.e.  they  only  require 
the  knowledge  of  the  tangent  stiffness  matrix.  However,  a  global  modal  analysis  based  on  the 
eigenvalues  of  the  full  system  is  computationally  prohibitive  as  the  problem  grows  larger.  This  fact 
makes  the  proposed  approach  attractive  for  the  study  of  the  global  instability  behavior  since  the 
computational  cost  is  of  the  same  magnitude  of  the  analytical  instability  criterion. 

In  the  next  section  we  summarize  the  research  accomplishments  of  this  project. 
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3  Research  accomplishments 

The  research  resulting  from  this  project  has  resulted  in  7  published  technical  paper  including  the 
following  major  contributions: 

•  A  formulation  that  can  capture  fracture  and  shear  banding  in  a  single  unified  model  [113] 

•  A  novel  hyperelastic  formulation  for  shear  bands  at  large  deformations  and  compared  it  with 
the  more  common  but  less  reliable  hypoelastic  formulation  [114] 

•  Stability  analysis  of  Shear  bands,  fracture  and  the  interaction  of  both  have  been  proposed  to 
provide  more  insight  into  the  point  of  instability  [1 15-1 17] 

•  The  unified  model  was  disceretized  using  finite  element  formulation  and  was  implemented 
into  a  research  finite  element  software  developed  in  the  group  [28] 

•  A  particular  ARL  application  of  interest,  the  failure  analysis  of  friction  stir  welded  (FSW) 
aluminum  joints  for  use  in  Army  land  vehicles,  has  been  studied  and  the  findings  were 
reported  [118] 

•  Numerical  simulations  on  one  and  two  dimensional  problems  show  improved  failure  model¬ 
ing  over  a  wider  range  of  strain  rates  than  either  fracture  or  shear  band  modeling  alone 

•  Benchmark  problems  (e.g.  plate  under  high  rate  tension,  the  Kalthoff  problem)  have  been 
investigated  and  a  brittle  to  ductile  transition  was  observed  which  is  in  line  with  theoretical 
and  experimental  observations. 

Next  we  summarize  in  more  details  some  of  the  research  accomplishments. 

3.1  A  unified  model  for  metal  failure  capturing  shear  banding  and  fracture 

3.1.1  Summary: 

In  this  work  a  thermodynamically  consistent  model  which  accounts  for  both  shear  banding  and 
dynamic  fracture  and  can  thus  capture  both  failure  bodes  at  intermediate  strain  rates,  is  presented. 
The  model  consists  of  an  elastic-viscoplastic  material  with  strain  hardening,  strain  rate  harden¬ 
ing,  and  thermal  softening.  Fracture  is  modeled  with  the  phase  field  method,  for  which  a  novel 
modification  is  presented  here  to  account  for  the  creation  of  fracture  surfaces  by  inelastic  work. 
Numerical  examples  are  presented  to  illustrate  the  basic  behavior  of  the  model,  and  to  compare  it 
to  three  special  cases:  a  damage  free  case,  an  isothermal  case,  and  an  isothermal  case  where  the 
contribution  of  inelastic  work  to  fracture  is  excluded. 

3.1.2  Dynamic  Fracture  Governing  Equations 

The  full  set  of  coupled  shear  band  and  fracture  PDFs  that  describe  dynamic  fracture  of  metals 
(ranging  from  low  to  high  impact  rates),  has  recently  been  proposed  by  Pl-Waisman  [113,  119], 
and  is  summarized  under  one  framework,  as  follows: 
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Conservation  of  Momentum:  poili  =  F^jTij^A  (1) 

Conservation  of  Energy:  PqcT  =  KJF^jF^jT^AB  +  X^9  (^)  T,  (2) 

2  /n  dTri/ 

Phase  Field  (damage)  Evolution:  c  =  4/qC  //  —  (fF’*'  +  P"*")  (3) 

Stress  Degradation:  Tjj  =  Tj”  +  m  (c)  (4) 

2 

Flow  Potential:  0  =  ^  —  1  =  0  (5) 

Inelastic  Constitutive  Relation:  ^  =  g{a,  T,  ?)  (6) 


Stress  Update:  £^bl  =  -2  6^;  (8) 


with  the  appropriate  boundary  and  initial  conditions.  Here  the  unknown  fields  are:  Ui  the 
displacements,  Tij  =  Joij  the  Kirchoff  stress  (cr  is  the  cauchy  stress  and  cr'  is  the  deviatoric 
component  of  the  cauchy  stress),  T  the  temperature,  c  the  phase  field  (damage)  parameter  and  P  the 
equivalent  plastic  strain.  Large  deformation  kinematics  is  accounted  for  through  the  deformation 
gradient  tensor  Fia  and  its  determinant  J,  where  the  initial  and  current  configurations  are  denoted 
by  upper  and  lower  case  indices,  respectively.  Note  that  the  formulation  builds  upon  a  hyperelastic 
formulation  (as  opposed  to  the  common  hypoelastic  modeling  approach)  as  presented  in  McAuliffe 
and  Waisman  [120]  and  Tij  is  derived  from  an  elastic  potential.  Eq.  (8)  is  an  incremental  kinematic 
update,  close  to  the  one  presented  by  Simo  and  Hughes  [121],  where  is  the  elastic  left  Cauchy 
Green  tensor,  and  dl,^  are  the  elastic,  plastic  and  thermal  rate  of  deformations,  and  £y  stands 

for  a  Lie  derivative  operator.  In  addition,  m(c)  is  the  phase  field  degradation  function,  1F+  and 
are  the  portion  of  the  elastic  and  plastic  strain  energies  being  degraded  by  fracture. 

The  particular  form  of  g  (a,  T,  P')  used  in  the  preliminary  work  is  a  thermo-visco-plastic  phe¬ 
nomenological  material  model  (a  modified  Litonski  model  [122]),  which  has  the  form 


Plastic  Deformation  Tensor: 


= 


ae 


d± 

dcr 


g{a,T,^) 


7o 


a 

CTo  (1  ^/7o)”  |l  -  5 

[exp(’’f*')  1] 

l}. 

(9) 


and  describe  the  strong  coupling  between  strain  hardening  and  temperature  softening.  More  details 
on  quantities,  material  parameters  and  discussion  on  the  physics  (not  presented  here  due  to  space 
limitations)  can  be  found  in  [123].  Shear-banding  is  captured  by  the  material  model  in  Eq.  (9) 
and  temperature  conductivity  (k  7^  0  in  Eq.  2)  is  used  to  introduce  a  weak  length  scale  into  the 
system  and  regularize  the  equations.  The  work  on  shear  bands  regularization  (without  cracks),  has 
been  funded  by  DOE,  where  it  was  proven  that  thermal  conductivity  leads  to  mesh  independent 
results,  see  for  example  [124-126].  In  the  current  proposal,  brittle  and  ductile  fracture  are  captured 
by  incorporation  of  the  phase  field  method,  as  described  in  Eq.  (3).  In  this  approach,  cracks  are 
described  as  diffusive  damage  entities  [127-129]  (closely  related  to  regularized  gradient  damage 
methods  [130]).  The  complete  multiphysics  system  accounts  for  plasticity,  heat  transfer,  brittle 
and  ductile  fracture  and  large  deformation  kinematics.  In  this  proposal  our  objective  is  to  enhance 
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the  ductile  fracture  component  by  adding  a  second  phase  field  equation  with  a  ductile  length  scale, 
which  is  driven  by  a  Gurson-Tvergaard-Needleman  (GTN)  void  growth  model  [131,  132]. 

Mixed  finite  element  formulation  has  been  proposed,  where  each  equation  is  multiplied  by  a 
corresponding  weight  function  and  integrated  over  the  domain.  This  weak  form  contains  spatial 
derivatives  of  displacement,  temperature  and  the  phase  field  parameter,  and  therefore  these  three 
fields  must  be  approximated  by  at  least  functions.  On  the  other  hand,  spatial  derivatives  of 
the  equivalent  plastic  strain  and  stresses  do  not  appear  and  thus  C~^  functions  will  suffice  (note 
that  a  Babuska-Brezzi  [133, 134]  condition  must  be  considered  when  developing  mixed  elements). 
The  weak  form  is  discretized  in  time  using  an  implicit  Newmark  method  where  at  each  time  step, 
a  nonlinear  problem  is  solved  until  convergence  for  the  five  fields  p  =  [uicT  aij  ,  using 
a  monolithic  Newton  scheme.  Such  numerical  approach  doesn’t  require  splitting  or  lagging,  as 
opposed  to  most  mesh  sensitive  explicit  schemes  used  for  these  problems.  Superior  behaviour  of 
the  implicit  scheme  in  terms  of  larger  time  steps,  stability  and  accuracy  has  been  reported  in  [135]. 
Consistent  linearization  of  the  nonlinear  system  leads  to  the  following  Jacobian  matrix  [1 19]: 


0 

0 

0 

Ma/3  +  +  C^a/3 

e'er 

'^a/3 

CCT 

0 

+  G^J 

C^Tcr 

^a(3 

+  La/3 

^cxl3 

r^aT 

'^a/3  +  S_a/3  + 

^a(3 

0 

0 

r^lpT 

M%  +  g; 

(10) 

where  M  denotes  mass  matrices,  K  stiffness  matrices  arising  from  linear  material  behavior  such  as 
the  Laplacian  terms  in  the  energy  and  phase  field  equations,  G  matrices  arising  from  the  derivatives 
of  the  flow  law  g  {a,  T,^).  L  matrices  arising  from  nonlinear  geometric  behavior  which  arise  due 
to  the  dependence  of  the  deformation  gradient,  its  inverse,  and  its  determinant  on  the  displacement, 
C  arise  from  derivatives  of  the  degradation  function  m  (c).  W  matrices  are  due  to  linearizing  W~^ 
with  respect  to  displacement.  Lastly,  H  is  due  to  linearizing  the  hyperelastic  material  law  defining 
the  stress.  The  structure  of  the  Jacobian  reflects  the  strongly  coupled  nature  of  the  system  of  PDEs, 
and  does  not  reveal  any  clear  way  to  simplify  or  reduce  the  system  by  eliminating  any  of  the  five 
fields  in  p.  Also  note  that  the  Jacobian  is  not  symmetric,  which  requires  appropriate  solvers. 

To  illustrate  the  potential  of  the  proposed  model  we  study  a  benchmark  problem  of  an  impact 
onto  a  notched  steel  plate  at  different  velocities,  also  known  as  the  Kalthoff  problem  [136,  137]. 
In  these  experiments,  which  are  somewhat  counterintuitive,  it  was  observed  that  above  a  critical 
impact  velocity,  the  failure  is  ductile,  with  a  shear  band  propagating  straight  from  the  notch  tip  in 
a  curved  fashion.  Below  the  critical  velocity,  a  crack  is  initiated  and  propagates  from  the  notch  tip 
at  about  70°.  Furthermore,  at  intermediate  rates,  a  brittle-to-ductile  failure  transition  is  observed 
where  both  failure  modes  are  present.  These  are  two  different  failure  phenomenon  since  shear 
bands  are  driven  by  shear  heating  caused  by  inelastic  deformations  while  cracks  are  driven  by 
tensile  loading.  We  emphasize  that  there  are  currently  no  models  that  can  capture  both  failure 
modes  at  the  same  time  with  a  single  model. 

Preliminary  results,  using  the  dynamic  fracture  model  proposed  by  the  PI,  are  shown  in  Figure 
5.  Projectile  impacts  at  10,  15,  20  and  25  m/s  are  presented.  For  each  impact  velocity,  we  show 
the  phase  field  parameter  and  the  equivalent  plastic  strain  pairs.  A  localized  phase  field  value 
close  to  1  indicate  a  crack  (brittle  fracture)  while  a  large  localization  of  plastic  strain  indicate  a 
shear  band  (ductile  failure).  We  consider  a  small  region  around  the  notch  tip  in  the  analysis  due 
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Figure  4:  Impact  experiments  on  a  notched  steel  plate  (Kalthoff  problem  [136,  137]).  (Left)  Low  velocity 
impact  leading  to  brittle  fracture  propagation  from  the  crack  tip  at  about  70°.  (Middle)  High  velocity  impact 
with  shear  band  forming  out  of  the  notch  tip  and  curving  inwards.  (Right)  Intermediate  impact  velocity 
leading  to  so  called  brittle-to-ductile  failure  transition  with  fracture  initiating  on  top  of  the  shearband. 

to  current  computational  limitations  on  the  size  of  the  problem.  Nonetheless,  the  model  is  able  to 
capture  both  shear  bands  and  cracks  at  different  strain  rates.  Furthermore,  it  is  interesting  to  note 
that  brittle  fracture  (at  low  impact  speed)  transitions  to  ductile  fracture  (at  higher  impact  speed) 
with  increase  in  strain  rate.  However,  at  this  point  the  analysis  is  terminated  as  the  nonlinear  solver 
diverges  due  to  severe  element  distortions  due  to  the  large  geometric  nonlinearities. 

The  aforementioned  formulation  requires  a  rather  fine  mesh  in  the  shear  band  or  crack  path  to 
be  effective.  In  general,  the  crack  and  shear  band  paths  are  not  known  a  priori,  and  therefore  an 
adaptive  method  is  needed  to  keep  computational  costs  reasonable  for  many  practical  situations. 
For  the  preliminary  simulations  here,  we  fully  refine  the  area  around  the  notch  tip,  but  can  only 
capture  the  initiation  of  failure. 

While  initial  promising  results  that  show  mesh  insensitivity  have  been  reported,  validation 
with  experiments  and  further  discoveries  are  limited  since  more  accurate  micromechanics  based 
constitutive  models  are  needed  for  ductile  fracture  modeling.  The  reader  may  refer  to  []  for  more 
details  on  model  derivation. 

3.2  Stability  analysis  of  Shear  bands  and  Fracture  at  high  strain  rate 

In  this  section,  we  analyze  the  physical  stability  of  both  the  fracture  and  the  shear  band  failure 
modes  and  their  interaction  using  a  linear  perturbation  method.  The  analysis  provides  insight  into 
the  dominant  failure  mode  and  can  be  used  as  a  criterion  for  mesh  refinement. 

Several  numerical  results  with  different  geometries  and  a  range  of  strain  rate  loadings  demon¬ 
strate  that  the  stability  criterion  predicts  well  the  onset  of  failure  instability  in  dynamic  fracture 
applications. 

For  the  example  problems  considered,  if  a  fracture  instability  precedes  shearbanding,  a  brittle- 
like  failure  mode  is  observed,  while  if  a  shear  band  instability  is  initiated  significantly  before 
fracture,  a  ductile-like  failure  mode  is  expected.  In  any  case,  fracture  instability  is  stronger  than  a 
shear  band  instability  and  if  initiated  will  dominate  the  response. 
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Figure  5:  Numerical  results  of  the  Kalthoff  experiment  in  Fig.  4  with  projectile  impacts  of  10, 
15,  20  and  25  mjs.  For  each  velocity  the  phase  field  parameter  and  the  equivalent  plastic  strain 
pairs  are  shown.  Note  how  brittle  fracture  (at  low  speeds)  transitions  to  ductile  fracture  (at  higher 
speeds)  with  increase  in  strain  rate. 

3.2.1  Problem  Statement 

The  dynamic  fracture  problem  studied  in  this  chapter  consists  of  a  thermo-mechanical  system 
in  which  fracture  is  modeled  through  the  phase-filed  formulation,  as  given  in  equation  set  (1)- 
(8).  Herein  we  assume  a  small-strain  formulation  for  elastic-viscoplastic  material  models  and 
in  addition  consider  the  ID  reduction  of  the  model.  Thus  the  strong-form  governing  equations 
degenerated  in  a  ID  pure  shear  formulation  from  the  multidimensional  set,  is  as  follows: 


Momentum: 

pii  =  r' 

(11) 

Damaged  Elastic  Const.  Law: 

dW^ 

t^GY  +  [m(c)  1] 

(12) 

Energy  balance: 

pCj,f  =  XT"  +  xrY 

(13) 

Strain-Displacement: 

—  u'  —  7^ 

(14) 

Phase-Eield: 

pe,c  =  GM'  -  -  ^(1L+  +  P+) 

z/q  UC 

(15) 

Inelastic  Const.  Law: 

T  =  PiT)Qi-,’’)R(Y) 

(16) 

where  G  is  the  shear  modulus  and  u'  is  the  total  strain.  Note  that  a  super-imposed  dot  {x)  corre¬ 
sponds  to  a  temporal  derivative  and  a  prime  {x')  corresponds  to  a  spatial  derivative. 

The  plastic  free  energy  that  contributes  to  fracture,  restricted  to  ID,  is  given  by 

rto 

P~^  =  /  (17) 

Jo 

Recalling  the  assumption  of  a  monotonic  loading  and  additionally  assuming  that  all  the  elastic 
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strain  is  contributing  to  W~^,  we  can  then  specify 


1  8W~^ 

W+  =  W  =  -G(rf  ^  =  GY  (18) 

Therefore,  the  damaged  elastie  constitutive  law  (12)  becomes 

T  =  m{c)GY  (19) 

Additionally  we  neglect  micro-inertia  effects  in  the  phase  field  equation  (15),  i.e.  6c  =  0,  and 
express  it  as 

c  =  ac"  -  (fT+  +  P+)  (20) 

with  a  =  4/q  and  /)  =  ^.  The  parameter  a  is  referred  to  as  the  gradient  coefficient  and  corre¬ 
sponds  to  the  square  of  the  characteristie  length  of  the  model.  As  mentioned  before,  2/o  defines  the 
width  of  the  diffused  crack,  which  means  that  as  /q  decreases,  the  value  of  a  also  decreases,  making 
the  crack  narrower.  The  parameter  (3,  relates  the  amount  of  energy  that  contributes  to  fracture  with 
the  critical  fracture  energy  Gc  (a  material  property).  Therefore,  the  term  /3(kF+  -I-  P+)  quantifies 
the  amount  of  energy  with  respect  to  Gc  driving  the  evolution  of  the  smeared  crack  and  serves  as 
the  souree  term  in  the  phase-field  equation. 

Applying  the  linear  perturbation  methodology  to  obtain  the  perturbed  equations,  we  get 


Momentum: 

Damaged  Elastie  Const.  Law: 

Energy  balance: 
Strain-Displaeement: 

Phase-Field: 

Inelastie  Const.  Law: 


p5u  =  5t'  (21) 

5t  =  G  -f  mo5Y^  (22) 

pc/T  =  XSr'  +  xStY  +  xroSY  (23) 

SY  =  Su'  -  5Y  (24) 


(ic  =  a5Y  -  /3^^W+{1  +  fl)5c  -  +  (iP+) 

(25) 


5t  =  -PoST  +  QoSY  +  RYY 


(26) 


where 


(27) 


5W+ 

5P+ 


dW^ 

'ryC 


dY 


SY 


(28) 

(29) 
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and 


Po  =  - 


dr 


X=Xo 


Qo  — 
Ro  = 


dr 

d^p 

dr 

d^p 


X=X0 


X=X0 


(30) 

(31) 

(32) 


Recall  that  a  variable  with  a  subscript  0  corresponds  to  the  value  of  that  variable  (or  its  deriva¬ 
tives)  at  the  current  equilibrium  point  which  is  being  perturbed.  For  example,  cq  corresponds  to 
the  value  of  the  phase-field  c  at  the  current  equilibrium  point.  The  spatial  and  time  derivatives  of 
the  independent  variables  correspond  to  the  multiplication  of  the  perturbation  of  the  independent 
variables  by  the  coefficients  ik  and  u,  respectively,  i.e.  Sx'  =  ik^a;  and  5x  =  coSx. 

Resolving  all  spatial  and  temporal  differentiations  through  the  aforementioned  procedure,  the 
independent  variables  can  be  eliminated  by  manipulation  of  (21)-(26)  and  the  characteristic  equa¬ 
tion  is  obtained.  Given  the  complexity  of  this  result  we  present  the  characteristic  equation  in  terms 
of  normalized  parameters,  following  a  similar  non-dimensionalization  as  suggested  by  Bai  [3]  for 
the  case  of  adiabatic  shear  bands,  that  is 


Cq  C\Cj  -\-  C/^Cj'^  —  0 


(33) 


with  the  following  expressions  for  the  coefficients  Go,  Gi,  (3*2,  G3  and  G4 
Go  =  Afk^C 

Ci  =  Afk^C  (l-5,  +  ^p) 

G2  =  Cg  [Af(  +  P^(l  —  0]  -|-  [Af(  (1  -|-  As)  -l- 1  -l-  P®(1  —  C)]  k^ 

G3  =  ^1  —  Ps  +  Agk^'j  +  [AyC  +  P®(1  —  C)] 

G4  =  ^ 
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where 


u  = 


A 

CpQo 


P  = 


pcjQo 


XToPp 

pCpQo 

X^Ppil 

pclQo 


(34) 

(35) 

(36) 

(37) 

(38) 


e- 


C  = 


a  I  (1+/P)  d’^my, 
^  ^  2  9c2  ' 


Af  = 


Gmn  Gmt) 


pa  = 


Qo 

Gmo  dP^ 


dr 

d^p 


'dW^ 


Qo  d-fP  /  dY 


(39) 

(40) 

(41) 


The  non-dimensional  regularization  parameter  (0),  which  accounts  for  the  effect  of  the  pertur¬ 
bation  length  in  relation  with  the  diffused  crack  width  (as  depicted  in  Figure  2),  is  defined  by 


0  =  1  +  =  1  +  {2lokf 


(42) 


and  the  non-dimensional  source  parameter  (0)  in  the  phase  field  equation  is  given  by 

=  2/3W+  =  2^VFo+  (43) 

\J  c 

Note  that  all  normalized  quantities  are  positive  with  the  exception  of  (^. 

To  study  the  characteristic  equation  in  (33)  we  employ  the  Routh-Hurwitz  stability  conditions[138, 
139],  which  state  that  a  fourth  order  polynomial  will  be  stable,  i.e.  all  its  root  are  in  the  negative 
real  half-plane,  if  the  following  conditions  are  met: 

Gi>Q  z  =  0,...,4  ,  CsCs  >  QCi  and  C'3C'2C'i  >  +  C'fC'o  (44) 


Therefore,  for  instability  to  occur  it  is  a  necessary  and  sufficient  condition  that  at  least  one  of 
these  inequalities  will  not  be  satisfied.  It  is  important  to  note  that  once  an  instability  condition 
is  met,  then  it  can  be  assumed  that  its  negation  is  true  for  all  subsequent  analysis.  Therefore, 
considering  the  limit  case  of  A:  — >  0,  there  are  two  conditions  for  which  the  characteristic  equation 
becomes  unstable,  namely 


where 


C<o 


GpF  —  0/0C  ~  1  >  0 


1  /  drriQ  \  ^  l+/o  d‘^mo 

mo  \  dc  )  2  dc^ 


(45) 

(46) 


15 


and 


CsB  —  1  —  Bg  <  0 


(47) 


where  Cpp  and  Csb  stands  for  the  instability  criterions  for  fracture  and  shear  bands,  respec¬ 
tively. 


Eigenvalue  Criterion:  The  stability  of  a  problem  can  also  be  analyzed  from  a  numerical  per¬ 
spective  by  making  use  of  the  concept  of  Lyapunov  stability  where  a  local  eigenvalue  analysis  of 
a  particular  partition  of  the  element  Jacobian  matrix  (or  Tangent  Stiffness  Matrix)  is  analyzed,  as 
presented  in  Arriaga  et  al.  [88]. 

As  Leroy  and  Ortiz  [140]  stated:  “arbitrarily  slow  perturbations  in  a  rate-dependent  solid  can 
only  grow  from  quasistatic  solutions.”  By  arbitrarily  slow  perturbations  it  is  understood  that  per¬ 
turbations  with  small  but  positive  growth-rate  or  the  first  eigenvalue  which  crosses  the  real  half 
plane  and  becomes  positive,  correspond  to  the  onset  of  instability. 

In  the  current  dynamic  fracture  model,  to  obtain  the  residual  equations  we  first  define  the 
weak  form  of  the  governing  equations  by  multiplying  each  equation  in  (l)-(8)  by  its  corresponding 
weight  function  [tw",  ,  integrating  over  the  problem  domain  and  using  integration 

by  parts  where  necessary. 


I’m  =  /  w'^pui  dO  -|-  /  w'^jaij  dO  —  /  w^aijTij  dT 

Jn  Jn  ’  Jr 


dn  -  / 

JQ  JQ 


ct‘S4,  +  Me)  - 1] 


ds^- 


yt  =  w^pCpT  dQ  +  /  w'^XT^i  dQ  —  /  7^)  dO  —  / uFXT^irii  dT 

Jn  Jn’  Jn  Jr 

Tc  =  j  w^c  dQ  +  j  dQ  +  j 

Jn  Jn’  Jn 


w'- 

n 


,2(o  dm 


(IL+  +  P+)  dO 


dfl-  w^^’giT,  d,  dO 

Jn  Jn 


(48) 

(49) 

(50) 

(51) 

(52) 


Note  that  the  Strain-Displacement  equation  is  not  explicitly  included  in  the  weak  form,  but  it  is 
implicitly  taken  into  account  when  computing  the  elastic  strain  in  the  Damaged  Elastic  Constitutive 
Law. 

Accounting  for  the  Babuska-Brezzi  condition  [141-143]  in  mixed  finite  element  formulations, 
the  shape  functions  for  each  field  must  be  chosen  with  care.  To  this  end,  we  choose  (7°  piecewise 
linear  shape  functions  for  displacement,  temperature  and  phase-field  parameter,  and  C~^  piecewise 
constant  functions  for  the  stress  and  equivalent  plastic  strain. 

The  residual  equations  can  be  grouped  into  a  residual  vector  r  and  a  solution  vector  x  with 
displacements,  temperature,  stress,  phase-field  and  plastic  strain  as 


Ui 

'r„‘ 

c 

rc 

T 

r  = 

Yt 

O'ij 

ra 

(53) 
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The  coupled  nonlinear  problem  can  then  be  stated  as 


— r(xo,  xo,  xo)  =  M  •  5x  +  C  •  (5x  +  K  •  (5x 


(54) 


where  r(x,  x,  x)  is  the  residual  vector  with  the  residual  of  each  equation  in  the  set  of  governing 
equations  of  the  problem  and  x  is  the  solution  vector  which  contains  all  field  variables  being  solved 
for  (Ui,  c,  T  Gij,  7^).  M-(ix,  C-(ix  and  K-(ix  are  the  obtained  by  computing  the  Gateaux  differential 
of  the  residual  (r)  in  the  (ix,  (ix  and  (ix  directions  as  follows 


M  •  (5x  =  d5xr(x,  X,  x) 
C  •  5x  =  d5xr(x,  X,  x) 
K  •  (ix  =  d5xr(x,  x,  x) 


— r(x,  X,  X  +  e5x) 
de 

-^r(x,  X  +  e(5x,  x) 
de 

^r(x  +  e(ix,  x,  x) 


€=0 


e=0 


e=0 


(55) 

(56) 

(57) 

(58) 


The  Jacobian  matrix  is  then  constructed  by  approximating  the  rate  and  acceleration  of  the  solu¬ 
tion  increment  (()x  and  ()x)  with  a  Newmark-beta  time-integration  scheme  such  that  (54)  becomes 

J  •  (ix  =  -r(xo,  Xo,  Xo)  (59) 

where  the  Jacobian  matrix  J  is  given  by 


'M* 

uu 

0 

0 

Kt, 

0 

C' 

^cu 

+  G,e 

GcT 

Gc(T 

Gcjp 

0 

0 

C^T  “1“  Gtt 

Gtct 

Gx^p 

(60) 

cru 

G(TC 

GctT 

Gq-^P 

0 

0 

GjPT 

GjPcr 

G^p;yP  +  GjPjP 

with  the  superscript  (*)  indicating  that  the  matrix  has  been  scaled  by  the  constants  that  result  from 
the  time-integration  scheme  used.  Here  K  =  -\-  G  represents  the  sum  of  the  stiffness  matrices 

associated  with  linear  material  behavior  such  as  elasticity  and  thermal  diffusion  (K^)  and  the 
tangent  stiffness  matrices  associated  with  material  non-linear  behavior  (G).  We  refer  to  McAuliffe 
and  Waisman  [27,  28]  for  additional  details  on  this  formulation. 

The  quasi-static  finite  element  formulation  of  the  system  can  be  acquired  by  considering  the  K 
component  of  the  Jacobain  matrix,  i.e.  the  component  of  the  Jacobian  matrix  that  is  affecting  the 
non-rate  terms  of  the  governing  equations.  Hence,  the  matrix  K  is  given  by 


■  0 

0 

0 

KL 

0 

^CU 

KJ;  +  Gee 

Gct 

Gccr 

Gq^P 

K  = 

0 

0 

H-XT  Gxt 

Gra 

Gx^p 

au 

G(tc 

GctT 

Gcr^P 

_  0 

0 

Gjpx 

GjPq- 

Gjp^p 

(61) 


Note  that  the  matrix  K  is  nonsymmetric  and  will  potentially  lead  to  complex  eigenvalues. 


17 


Therefore,  the  stability  condition  based  on  the  numerical  approximation  is  obtained  when  there 
exists  an  eigenvalue  of  K  with  a  positive  real  part,  i.e. 


Re  [eig  (K)]  >  0 


(62) 


The  main  advantage  of  this  process  is  that  it  is  straight  forward  to  apply  and  will  be  as  accurate 
as  the  FEM  approximation.  An  implicit  assumption  of  this  method  is  that  the  values  for  the  wave¬ 
length  of  the  perturbation  are  limited  to  permutations  of  the  degrees  of  freedom  within  the  element 
discretization.  In  other  words,  the  wave-lengths  are  locked  to  the  element  size.  Nonetheless,  we 
neglect  the  influence  of  this  effect  since,  as  mentioned  earlier,  the  limit  case  is  given  by  A:  — )•  0, 
which  corresponds  to  an  infinitely  large  wave-length  of  the  perturbation  (i.e.  uniform  solution), 
which  can  be  easily  captured  by  the  Finite  Element  discretization  of  the  problem. 

While  this  approach  will  be  shown  in  the  next  section  to  agree  well  with  both  analytical  criteria 
for  shear  bands  and  fracture,  one  limitation  of  this  method  is  that  it  does  not  distinguish  between 
the  types  of  failure  modes  as  is  clearly  obtained  with  the  analytical  criterion.  Another  drawback  of 
the  spectral  method  is  its  higher  computational  burden,  as  it  requires  a  solution  of  an  eigenvalue 
problem  at  each  element  and  will  slow  down  the  analysis,  which  might  be  significant  in  large  scale 
parallel  computations. 


3.2.2  Numerical  Results 

Square  in  tension:  In  this  section  we  study  the  analytical  criteria  on  a  2D  formulation  of  an 
initially  homogeneous  system,  where  the  effects  of  shear  band  localization  and  fracture  are  present. 
A  homogeneous  system  is  one  where  all  quantities  are  initially  constant  or  vary  smoothly.  This  has 
the  advantage  that  the  cause  for  localized  behavior  is  mostly  due  to  the  instability  of  the  material 
and  not  due  to  boundary  conditions  or  geometry. 

A  steel  square  plate  is  considered  as  shown  in  Figure  1  with  dimension  H=100//m.  The  plate 
is  stretched  uniaxially  under  a  displacement  control  loading  with  i)r=0-005m/s  and  5m/s,  which 
corresponds  to  a  nominal  strain-rate  of  5  x  10^  s“^  and  5  x  10"^  s“k 

A  hyperbolic  secant  type  imperfection  is  used  to  scale  the  yield  stress,  the  hardening  parameter 
and  the  critical  fracture  energy  by  a  constant 

The  mesh  consists  of  50  x  50  quad  elements,  modeled  by  an  irreducible  mixed-finite  element 
formulation  with  a  B-bar  implementation  to  reduce  shear-locking  effects.  As  defined  before,  the 
material  considered  is  a  4340  Steel  with  a  Johnson-Cook  constitutive  law  [144]  with  the  effect  of 

and  we  take  into  account  the  non-linear  effects  of  the  Taylor-Quinney  coefficient  as  proposed 
by  Vural  et  al.  [19],  which  takes  into  account  the  influence  of  the  isothermal  conditions  usually 
verified  in  quasi-static  loadings. 

Figure  6  shows  the  stress  strain  curve  of  the  two  different  strain-rates.  The  stress  is  the  average 
(722  in  the  entire  plate.  As  expected  from  the  Johnson-Cook  material  law,  a  higher  value  of  the 
strain-rate  increases  the  yield  stress  of  the  material.  Note  that  at  this  point  it  is  not  clear  what  type 
of  collapse  is  dominating  in  each  case. 

Because  a  homogeneous  problem  is  considered,  the  length-scales  are  large  in  comparison  with 
the  dimension  of  the  plate  being  represented.  This  causes  the  stability  conditions  to  be  less  local¬ 
ized  than  they  would  otherwise  be  in  a  problem  where  the  length-scales  are  considerably  smaller 
than  the  problem  size.  This  will  be  seen  in  the  Kalthoff  example  in  the  next  section. 
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Nominal  Strain 


Figure  6:  Stress-strain  curves  for  both  strain-rates.  The  low  strain-rate  corresponds  to  a  value  of 
5  X  10^  s“^  and  the  high  strain-rate  corresponds  to  a  value  of  5  x  10^  s“^. 

Low  Strain-rate  In  Figure  7  we  show  the  stress  strain  curve  for  the  low  strain-rate  as  well  as  a 
few  selected  snapshot  points  for  which  a  full  depiction  of  the  plate  will  be  shown.  The  repeated 
values  for  the  times  are  due  to  the  fact  that  the  failure  process  is  very  fast  and  the  difference  between 
two  points  in  time  is  extremely  small  and  approximately  0.2//S. 


Nominal  Strain 


Figure  7:  Stress-strain  curve  with  the  different  snapshot  points  represented. 

Figures  8-11  represent  the  time  snapshots  reporting  the  Phase-field  parameter  (c),  the  equiv- 
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alent  plastic  strain  (7^’)  and  the  instability  conditions  for  fracture  Cpp  (Eq.  45)  and  shear  bands 
CsB  (Eq.  47).  A  material  region  which  remains  stable  during  loading  is  marked  by  a  white  color 
while  a  colored  region  implies  instability  due  to  fracture  or  shear  banding.  Different  colors  show 
the  distribution  of  the  actual  value  of  the  instability  conditions  and  their  intensity. 

The  first  snapshot  in  Eigure  8  is  taken  near  the  peak  stress  and  it  can  be  seen  that  the  problem 
only  shows  an  instability  based  on  a  shear-band  criterion.  The  threshold  of  -1  for  Csb  delineated 
the  region  where  the  equivalent  plastic  strain  is  concentrating  and  starting  to  localize,  although  this 
effect  is  not  significant  as  the  plate  has  a  more  or  less  uniform  distribution  of  the  equivalent  plastic 
strain. 

The  last  three  snapshots  (Eigures  9-11)  show  the  onset  of  fracture,  and  its  propagation  until 
the  complete  rupture  of  the  solid.  It  is  interesting  to  note  that,  as  the  instability  condition  for 
fracture  precedes  and  delineates  the  crack,  the  shear  band  instability  disappears.  This  is  a  natural 
consequence  of  the  unloading  that  occurs  in  the  vicinity  of  the  crack.  It  is  also  interesting  to  note 
that  fracture  propagates  in  mode  I  along  the  boundary,  as  would  be  expected  at  this  slow  strain  rate. 

High  Strain- rate  In  Figure  12  we  show  the  stress  strain  curve  for  the  high  strain-rate  as  well  as 
a  few  selected  snapshot  points  for  which  a  full  depiction  of  the  plate  will  be  shown. 

Similarly  to  the  previous  section.  Figures  13-16  represent  the  Phase-field  (c),  the  equivalent 
plastic  strain  (7^)  and  the  instability  conditions  for  fracture  Cpp  (Eq.  45)  and  shear  bands  Csb 
(Eq.  47).  For  the  stability  conditions,  a  white  region  signifies  that  the  material  is  stable  while  a 
colored  region  implies  instability.  Different  colors  show  the  distribution  of  the  actual  value  of  the 
stability  conditions  and  their  intensity. 

The  first  snapshot  in  Figure  13  shows  a  situation  that  is  similar  to  the  initial  snapshot  of  the 
previous  section.  Here,  a  shear  band  instability  is  visible,  but  unlike  in  the  previous  example,  the 
shear  band  is  more  developed  and  greater  inhomogeneity  in  the  equivalent  plastic  strain  is  already 
visible.  The  threshold  of  -1  for  Csb  delineates  well  where  the  localization  of  plastic  deformation 
is  most  pronounced.  Fracture  is  virtually  non-existent  with  only  a  slight  homogeneous  increase  in 
the  phase  field  parameter  and  no  fracture  instability. 

The  second  snapshot  in  Figure  14  shows  a  much  more  developed  shear  band.  The  equivalent 
plastic  strain  is  localized  in  the  45*^  direction  and  the  Csb  stability  condition  demarcates  the  shear 
band  location.  A  Cpp  instability  is  already  visible  in  the  regions  with  large  plastic  deformation 
owing  to  the  effect  of  P^. 

In  the  third  snapshot  represented  in  Figure  15,  two  cracks  nucleate  in  the  corners  of  the  plate 
where  the  plastic  deformation,  and  consequently  the  value  of  P+,  are  the  largest.  It  is  interesting  to 
note  that  the  cracks  initially  propagate  in  mode  I  as  can  be  seen  in  the  distribution  of  the  phase-field 
and  in  the  fracture  stability  condition.  At  the  same  time,  the  shear  band  stability  condition  reverts 
back  into  a  stable  state  near  the  cracks,  since  this  region  has  unloaded,  but  the  center  of  the  plate 
is  still  uncracked  and  unstable  from  the  shear  band  point  of  view. 

Finally,  the  last  snapshot  is  shown  in  Figure  16,  where  the  crack  has  propagated  throughout  the 
plate  causing  the  complete  failure  of  the  solid.  The  fracture  stability  condition  recovers  the  full 
geometry  of  the  crack  and  the  shear  band  stability  condition  has  completely  subsided  back  into 
the  stable  regime.  It  is  also  interesting  to  note  that  fracture  appears  on  top  of  the  shear  band  and 
dominates  the  failure  of  the  plate. 
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Figure  8:  Low  Strain-rate,  t  =  22.0ms. 
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(c)  Fracture  instability 
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Figure  9:  Low  Strain-rate,  t  =  23.7ms  (A). 
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Figure  10:  Low  Strain-rate,  t  =  23.7ms  (B). 
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(c)  Fracture  instability 


CsB 

(d)  Shear-band  instability 


Figure  11:  Low  Strain-rate,  t  =  23.7ms  (C). 
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Figure  12:  Stress-strain  curves  with  different  snapshots.  The  repeated  values  for  the  times  are  due 
to  the  fact  that  the  storage  of  the  time  does  not  have  enough  precision  to  capture  the  increment 
between  each  step  of  0.2//s. 


Kalthoff  problem  In  this  section  we  study  the  analytical  criteria  on  a  2D  formulation  of  the 
Kalthoff  problem[145]  with  the  geometry  given  in  Figure  17. 

The  material  is  a  representative  steel  modeled  through  a  combined  modified  version  of  the 
models  of  Zhou  et  al.  [21],  Areias  and  Belytschko  [146,  147]  and  described  in  McAuliffe  and 
Waisman  [28].  The  flow  law  is  therefore  given  by 


T,  Y)  =  To 


a 

(Jo(l  +  7^/7^)” 


(63) 


where 

with 


r(T)  =  max  [l  —  (i(e®  —  1),  1  —  tanh((iaft0 


T  —  To 

0  =  — - —  and  Sab  = 


ln(2) 


(64) 

(65) 


The  material  parameters  used  are  given  in  Table  1,  with  a  value  of  Gc  =  2.2E4  kJ/m^. 

Three  different  velocities  (uq)  are  applied  to  the  plate  (15,  20  and  25  m/s)  and  the  velocity  is 
applied  gradually  over  a  period  of  0.5/us.  Henceforth  these  three  velocities  will  be  referred  to  as 
Low,  Intermediate  and  High  velocity,  respectively. 

The  results  are  shown  only  in  the  vicinity  of  the  notch  tip,  where  the  mesh  was  significantly 
refined.  Specifically,  the  figures  represent  an  area  of  1  x  Imm^  and  the  lower-left  corner  of  the 
figures  will  be  at  x=50.4mm  and  y=l  12.9  mm. 
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Figure  13:  High  Strain-rate,  t  =  14.6/iS. 
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Figure  14:  High  Strain-rate,  t  =  20.0/iS. 
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Figure  15:  High  Strain-rate,  t  =  22.4/is. 


28 


Figure  16:  High  Strain-rate,  t  =  23.8/iS. 
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THICKNESS:  1/4”  (6.35mm) 


Figure  17:  Full  geometry  of  the  Kalthoff  Problem. 


Table  1 :  Areias-Belytsehko  parameters 


Property  Name 

Symbol 

Value 

Unit 

Reference  strain-rate 

Yo 

1.0 

10-3/s 

Strain-rate  hardening  parameter 

m 

70 

- 

Yield  stress 

(^0 

2000.0 

MPa 

Yield  strain 

7o 

0.01 

- 

Strain  hardening  exponent 

n 

0.01 

- 

Reference  Temperature 

To 

293 

K 

Thermal  softening  parameter 

S 

0.8 

- 

Thermal  softening  parameter 

k 

500 

K 
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Low  Impact  Velocity:  Two  different  times  are  represented  for  this  case:  27.40//S  and  27.84;us. 
Figures  18-19  report  the  Phase-field  parameter  (c),  the  equivalent  plastic  strain  (7^)  and  the  in¬ 
stability  conditions  for  fracture  Cpp  (Eq.  45)  and  shear  bands  Csb  (Eq.  47).  For  the  stability 
conditions,  a  white  region  signifies  that  the  material  is  stable  while  a  colored  region  implies  insta¬ 
bility.  Different  colors  show  the  distribution  of  the  actual  value  of  the  stability  conditions  and  their 
intensity. 

The  first  snapshot  in  Figure  18  corresponds  to  a  point  immediately  before  the  crack  starts 
propagating.  Both  the  fracture  instability  and  the  shear  band  instability  are  clearly  visible  as  both 
mechanisms  of  failure  start  competing.  Some  plastic  deformation  is  visible  in  the  bottom  of  the 
notch  where  the  shear  band  wants  to  develop  and  the  value  of  the  phase  field  is  clearly  nucleating 
a  crack  on  the  upper  part  of  the  notch. 

The  next  snapshot  shown  in  Figure  19  illustrates  an  advanced  stage  of  crack  propagation  .  The 
fracture  instability  identifies  the  crack  and  its  propagation  direction  and  the  shear  band  condition 
recovers  stability  due  to  the  unloading  of  the  crack.  At  this  impact  velocity,  the  contribution  of 
thermal  diffusion  stabilizes  the  onset  of  the  shear  band  and  this  mechanism  is  unable  to  fully 
develop  before  the  crack  nucleation  and  growth. 

Intermediate  Impact  Velocity:  Two  different  times  are  represented  for  this  case:  25.21;us  and 
25.66/iS.  Figures  20-21  represent  for  this  problem  the  same  quantities  as  before. 

The  first  snapshot  in  Figure  20  corresponds  to  a  point  immediately  before  the  crack  starts 
propagating.  Similarly  to  the  first  snapshot  in  the  Low  impact  velocity  case,  the  fracture  instability 
is  clearly  visible  as  well  as  the  shear  band  instability.  However,  the  shear  band  instability  condition 
is  now  much  more  pronounced,  with  values  going  below  the  -1.5  threshold,  and  the  shear  band  is 
clearly  visible  in  the  plot  of  the  equivalent  plastic  strain.  At  this  point  the  shear  band  localization 
is  underway  and  fracture  is  about  to  initiate. 

The  last  snapshot  in  Figure  21  shows  an  advanced  stage  of  the  propagation  of  the  crack.  The 
fracture  instability  identifies  the  crack  and  its  propagation  direction  and  the  shear  band  condition 
has  now  significantly  receded  due  to  the  unloading  of  the  crack.  The  appearance  of  the  crack  has 
clearly  halted  the  continuing  development  of  the  shear  band  as  can  be  seen  in  the  value  of  the 
equivalent  plastic  strain,  which  has  not  changed  noticeably  since  the  previous  state.  Nonetheless, 
this  intermediate  strain  represent  a  case  in  which  both  fracture  and  shear  banding  processes  are 
present  at  the  same  time. 

High  Impact  Velocity  Two  different  times  are  represented  for  this  case:  19.66/us  and  21.75/us. 
Figures  22-23  represent  for  this  problem  the  same  quantities  as  before. 

The  first  snapshot  in  Figure  22  is  taken  at  a  moment  where  the  development  of  the  shear  band 
is  at  a  similar  stage  than  in  the  first  snapshot  of  the  Intermediate  velocity  (Fig.  20).  However,  due 
to  the  higher  rate  of  deformation,  the  shear  band  is  steadily  progressing  and  the  onset  of  fracture 
instability  is  barely  present.  At  this  point,  and  unlike  the  previous  velocity,  there  is  not  a  strong 
indication  that  a  crack  is  nucleating. 

The  final  snapshot  in  Figure  23  shows  an  advanced  stage  of  the  propagation  of  the  shear  band. 
The  shear  band  stability  condition  is  strongly  demarcating  the  advance  of  the  shear  band,  partic¬ 
ularly  at  the  -1.5  threshold  and  the  fracture  stability  condition  is  only  showing  mild  instability. 
There  is  no  apparent  sign  of  the  nucleation  of  a  crack  an  the  failure  mechanism  seems  to  be  fully 
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Figure  18:  Low  Velocity,  t  =  27.40/iS. 
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Figure  19:  Low  Velocity,  t  =  27.84/iS.  Elements  belonging  to  the  crack  (i.e.  c  >  0.97)  were 
removed. 
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(c)  Fracture  instability  (d)  Shear-band  instability 

Figure  20:  Intermediate  Velocity,  t  =  25.21/iS.  Elements  belonging  to  the  crack  (i.e.  c  >  0.97) 
were  removed. 
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Figure  21:  Intermediate  Velocity,  t  =  25.66/iS.  Elements  belonging  to  the  crack  (i.e.  c  >  0.97) 
were  removed. 
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dominated  by  shear  band  localization. 
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Figure  22:  High  Velocity,  t  =  19.66/iS.  Elements  belonging  to  the  crack  (i.e.  c  >  0.97)  were 
removed. 
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Figure  23:  High  Velocity,  t  =  21.75/iS.  Elements  belonging  to  the  crack  (i.e.  c  >  0.97)  were 
removed. 
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3.3  Application  the  failure  analysis  of  friction  stir  welded  (FSW)  aluminum 
joints  for  use  in  Army  land  vehicles 

3.3.1  Summary: 

Prediction  of  shear  band  formation  and  other  strain  localization  processes  present  many  computa¬ 
tional  challenges  that  must  be  overcome  to  enable  dynamic  failure  prediction  and  material  design 
of  ductile  material  systems.  The  current  work  presents  a  finite  element  based  computational  frame¬ 
work  accounting  for  this  critical  deformation  process,  as  applied  to  a  detailed  investigation  of  fric¬ 
tion  stir  welded  (FSW)  aluminum  joints.  A  stir  welded  joint  has  several  zones,  each  with  distinct 
microstructural  characteristics  and  material  properties.  For  applications  in  Army  land  vehicles, 
which  may  be  subject  to  under-body  blast,  an  understanding  of  the  energy  absorption  capability  of 
these  joints  is  needed.  Thus  material  inhomogeneity,  dynamic  loading,  and  detailed  understand¬ 
ing  of  small  scale  failure  processes  must  all  be  accounted  for  to  accurately  model  FSW  material 
behavior.  In  this  study,  an  implicit  nonlinear  consistent  (INC)  or  monolithic  solution  technique  is 
used  to  predict  shear  band  formation  and  estimate  the  energy  absorption  and  failure  strain  of  a  stir 
welded  aluminum  joint.  It  has  been  shown  that  failure  initiating  at  material  interface  regions  can 
be  predicted,  and  furthermore  that  abrupt  material  property  gradients  predominantly  contributes  to 
FSW  joint  failure. 

3.3.2  Friction  Stir  Welding  (FSW) 

High  strength  lightweight  aluminum  alloys  offer  potential  advantages  as  replacements  for  tradi¬ 
tional  steels  in  many  Army  vehicles  in  terms  of  weight  specific  mechanical  properties.  Unibody 
chassis  construction,  as  opposed  to  body  on  frame  construction,  is  being  pursued  to  lend  enhanced 
rigidity  to  maintain  structural  integrity  during  potential  under-body  blast  events.  Since  unibody 
construction  requires  the  elimination  of  bolted  joints,  weldability  of  the  chassis  material  is  cru¬ 
cial.  Alloys  from  the  aluminum  2XXX  and  5XXX  series  are  notoriously  difficult  to  weld  with 
conventional  techniques,  but  can  be  joined  with  Friction  Stir  Welding  (FSW)  [148]. 

FSW  is  a  solid  state  joining  process,  which  has  recently  been  shown  to  be  capable  of  producing 
joints  in  aluminum  up  to  76.2  mm  thick  [149,  150].  The  FSW  tool  consists  mainly  of  a  shank, 
shoulder,  and  pin  (shown  in  Figure  24),  which  rotate  as  they  advance,  stirring  the  material  together. 
Significant  inelastic  deformation,  heat  production,  and  dynamic  recrystallization  occur  during  this 
process,  which  results  in  the  formation  of  several  zones  with  distinct  microstructure  and  material 
properties  [151-156]. 

Figure  25  shows  a  cross-section  of  a  typical  FSW  joint,  in  this  case  2139  Aluminum,  which 
illustrates  the  nature  of  these  distinct  material  zones.  The  material  zone  most  immediately  sur¬ 
rounding  the  high-torque  tool  pin,  inserted  between  the  welded  plates,  is  characterized  by  an  upper 
and  lower  weld  nugget.  The  upper  and  lower  weld  nuggets  are  zones  B  and  A  respectively  in 
Figure  25.  This  comes  both  as  a  direct  influence  of  the  tool  pin  on  the  material  grain  structure  as 
well  as  the  thermal  properties  of  the  material  and  recrystallization  processes.  The  tool  shoulder 
creates  a  great  deal  of  friction  in  contact  with  the  plate  surface  during  the  FSW  process  which  gen¬ 
erates  an  inordinate  amount  of  heat  in  comparison  to  the  bottom  of  the  plate  [157].  The  difference 
in  thermal  inputs  and  boundary  conditions  thus  creates  a  through  thickness  variation  in  recrys¬ 
tallization  and  grain  growth  which  leads  to  the  distinct  lower  and  upper  weld  nugget,  marked  as 
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zones  A  and  B,  respectively  in  Fig.  25.  Note,  that  for  thin  welds  there  is  no  readily  discernible 
difference  between  the  upper  and  lower  nuggets,  however  welds  of  at  least  2  cm  thick  are  of  the 
greatest  interest  for  Army  land  vehicle  applications.  At  this  thickness  there  is  a  clear  distinction 
between  the  upper  and  lower  nuggets,  as  shown  in[158,  159]  Just  outside  of  the  weld  nugget,  zone 
C  in  Figure  25,  is  a  third  material  zone  which  is  still  subject  to  both  the  mechanical  and  thermal 
influence  of  the  FSW  stirring  pin.  Further  from  the  tool  is  a  fourth  material  zone,  too  distant  to 
be  influenced  by  mechanical  stirring,  but  still  subject  to  thermal  microstructural  effects  as  heat  is 
conducted  away  from  the  FSW  joint  during  processing.  This  is  known  as  the  thermal  affected  zone 
and  is  marked  by  zone  D  in  Figure  25.  The  importance  of  accounting  for  these  variations  must 
be  noted,  and  they  will  play  a  significant  role  in  the  simulations  below.  Computational  modeling 
of  the  FSW  process,  aimed  at  numerically  predicting  the  FSW  joint  zones,  have  been  carried  out 
by  [160-162]  on  joints  6.35  mm  thick  and  by  [163-165]  for  joint  7  millimeters  thick.  Additional 
computational  studies  by  [166,  167]  sought  to  determine  how  the  FSW  process  parameters  such 
as  tool  advancement  speed  affect  the  final  weld.  The  studies  cited  above  used  either  a  Lagrangian 
or  Arbitrary  Lagrangian  Eulerian  (ALE)  type  of  numerical  formulation  with  adaptive  remeshing. 
However,  Eulerian  formulations  are  also  available  [168]. 

In  this  study,  we  numerically  model  high  rate  loading  of  a  small  cross  section  of  a  stir  welded 
joint.  The  goal  is  to  develop  a  predictive  capacity  for  energy  absorption  and  failure  of  FSW  joints 
in  dynamic  loading,  as  well  as  to  gain  an  understanding  of  the  relationship  between  an  FSW  joint 
microstructure  and  the  resulting  energy  absorption  capability  and  failure  strain.  It  is  found  that  the 
presence  of  the  abruptly  changing  material  zones  leads  to  significantly  reduced  energy  absorption 
capacity  of  the  FSW  joint  in  comparison  to  two  benchmark  tests  with  uniform  material  properties. 
The  benchmarks  used  the  material  properties  of  the  unaffected  material  and  the  lower  nugget,  the 
nominally  weakest  of  the  weld  zones.  This  suggests  that  weld  capacity  can  be  improved  by  altering 
the  FSW  process  parameters  such  that  the  final  weld  contains  more  smoothly  varying  zones 

3.3.3  Representative  Results 

Cross  weld  tension  Results  of  the  cross  weld  tension  simulations  at  two  different  strain  rates  are 
shown  below  in  Figures  26-27.  Figure  26  shows  characteristic  contour  plots  of  equivalent  plastic 
strain  at  the  time  of  failure,  and  Von-Mises  stress  at  the  time  of  peak  VV.  Contour  plots  are  shown 
side  by  side  for  a  uniform  region  of  untreated  aluminum  followed  by  a  plot  of  the  simulated  weld 
zone  region  in  Figure  26.  The  shear  band  pattern  for  the  UA  and  SW  cases  are  clearly  different, 
with  the  shear  band  in  the  SW  case  being  contained  entirely  in  zones  1  and  2.  For  the  stir  weld  case. 
Yielding  is  first  observed  in  the  zone  2,  which  corresponds  to  the  lower  weld  nugget.  Localization 
is  also  first  observed  in  this  zone,  but  shear  band  initiation,  as  indicated  by  a  precipitous  drop 
in  the  local  stress  level  during  loading,  originates  from  the  upper  weld  nugget.  This  behavior  is 
observed  for  both  loading  rates  in  the  numerical  simulation.  Quasi  static  experiments  on  FSW 
samples  reported  in  [158,  159]  showed  that  the  failure  origin  zone  would  be  determined  by  which 
zone  reaches  strain  hardening  saturation  first,  which  was  found  to  be  the  LN.  Herein  we  have 
matched  these  baseline  quasi-static  behaviors  and  extrapolated  both  to  higher  strain  rates  and  a 
greater  number  of  loading  conditions.  The  dynamic  case  is  somewhat  more  complicated,  with 
rate  effects  and  thermal  softening  being  present.  Strain  hardening  saturation  is  not  an  adequate 
indicator  for  failure  because  strain  softening  and  localization  does  not  necessarily  imply  failure, 
which  has  been  noted  in  the  experimental  work  of  Marchand  and  Duffy  [169],  and  is  evident  in 
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Figure  24:  FSW  tool.  The  tool  shoulder  is  butted  against  the  workpiece  and  the  tool  is  rotated  as  it 
advances.  This  process  mechanically  stirs  the  workpiece  metals  together  at  elevated  temperature, 
but  without  melting. 

the  numerical  results  presented  here.  The  nominal  strain  to  failure,  strain  energy  at  failure,  and  the 
zone  in  which  failure  initiates,  as  shown  in  Table  2. 

Perhaps  the  most  important  feature  of  Figure  26  is  the  buildup  of  stresses  that  occur  near 
the  material  zone  boundary.  This  is  essentially  an  interfacial  stress  directly  due  to  the  material 
inhomogeneity  and  property  mismatch,  which  leads  to  failure  of  the  SW  joint  at  lower  strains 
compared  to  the  unaffected  aluminum.  This  effect  leads  to  significant  differences  in  the  overall 
stress  strain  response  of  the  joint,  which  is  shown  in  Figure  27,  where  the  domain  averaged  stress 
(722  is  plotted  against  the  nominal  strain.  The  oscillations  present  in  these  plots  are  due  to  elastic 
wave  reflections,  which  are  damped  quickly  by  plastic  deformation.  This  illustrates  that  in  an 
averaged  sense,  the  FSW  joint  is  marginally  stronger  than  the  lower  weld  nugget,  but  at  the  cost  of 
significant  ductility,  failing  at  less  than  half  the  nominal  strain  of  the  untreated  aluminum.  The  loss 
of  mechanical  capacity  of  the  weld  can  clearly  be  seen  in  comparing  the  FSW  weld  performance 
(SW)  to  that  of  unwelded  aluminum  (UA)  in  Figure  28.  The  FSW  joint  fails  at  approximately 
5  times  lower  mechanical  energy  for  the  strain  rate  of  1E3  s~^  and  3  times  lower  for  the  rate  of 
1E4  due  both  to  softening  (and  lower  initial  rate  of  strain  energy)  and  a  lower  strain  capacity 
before  failure  initiation  see  Table  2.  Though  thermo-mechanical  property  degradation  due  to  FSW 
processing  clearly  weaken  the  FSW  joint  strength,  it  is  readily  apparent  that  inhomogeneity  of 
the  joint  and  interfacial  stress  risers  are  the  key  drivers  for  loss  of  mechanical  capacity.  This  is 
evidenced  by  the  fact  that  if  the  entire  region  were  comprised  of  the  nominally  weakest  material 
of  the  lower  weld  nugget,  this  would  still  outperform  the  actual  FSW.  This  comparison  is  clearly 
seen  in  Figure  28,  in  which  the  rate  of  strain  energy  is  comparable,  yet  the  duration  of  significant 
energy  absorption  prior  to  failure  is  much  larger  for  monolithic  LN  material  as  compared  to  the 
actual  FSW  joint.  The  reduction  in  energy  absorption  at  the  strain  rate  of  1E3  is  greater  than 
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Figure  25:  FSW  Cross  section  for  2139  Aluminum.  The  labeled  areas  correspond  to  different 
weld  zones  which  have  distinct  microstructural  and  mechanical  properties.  Zone  A  is  the  lower 
weld  nugget,  B  the  upper  weld  nugget,  C  the  thermo-mechanical  affected  zone,  and  D  the  thermal 
affected  zone. 

the  reduction  at  the  strain  rate  1E4  This  indicates  that  the  effect  of  material  inhomogeneity 
decreases  with  increasing  strain  rate.  This  is  due  to  the  fact  that  the  material  zones  have  the  same 
rate  hardening  characteristics,  which  become  increasingly  prevalent  relative  to  the  inhomogeneous 
strain  hardening  characteristics  as  the  applied  strain  rate  increases. 
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(a)  Unaffected  Aluminum  at  failure,  Nominal  Strain  Rate  (b)  Stir  Weld  at  failure,  Nominal  Strain  Rate  1E3 
1E3 


(c)  Unaffected  Aluminum  at  peak  VV,  Nominal  Strain  Rate  (d)  Stir  Weld  at  peak  VV,  Nominal  Strain  Rate  1E3 
1E3 

Figure  26: 
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(a)  Tensile  Stress  in  the  Direction  of  Loading,  Nomi-  (b)  Tensile  Stress  in  the  Direction  of  Loading,  Nomi¬ 
nal  Strain  Rate  1E3  nal  Strain  Rate  1E4 

Figure  27:  Domain  averaged  tensile  stress  in  the  direction  of  loading  vs  nominal  strain  for  each 
loading  case.  The  x  on  the  plot  indicates  the  failure  point. 


Nominal  Strain  Nominal  Strain 

(a)  W,  Nominal  Strain  Rate  1E3  (b)  VV,  Nominal  Strain  Rate  1E4 


Figure  28:  Rate  of  stress  working  for  each  case.  The  x  on  the  plot  indicate  the  failure  point 


Load  Config.  ^ 

Cross  Weld  Tension 

Nom  strain  rate  ^ 

1E3 

Mat  Config. 

Nom.  Failure  Strain 

W  at  Failure 

Failure  Origin 

UA 

0.48 

2.64e-i-08 

N/A 

LN 

0.75 

2.58e-l-08 

N/A 

SW 

0.15 

5.23e-l-07 

UN-TMAZ  Interface 

Nom  strain  rate  — )■ 

1E4 

Mat  Config.  1 

Nom.  Failure  Strain 

W  at  Failure 

Failure  Origin 

UA 

0.41 

2.24e-l-08 

N/A 

LN 

0.75 

2.57e+08 

N/A 

SW 

0.22 

7.62e+07 

UN-TMAZ  Interface 

Table  2:  Results  summary  for  the  cross  weld  tension  test:  nominal  strain  at  failure,  W  at  failure, 
and  failure  initiation  zone  ae  tabulated 
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4  Research  Dissemination 

Research  Results  were  disseminated  through  conference  presentations,  invited  seminars  and  keynote 

lectures,  posters  and  refereed  journal  publications. 

Resulting  Journal  Publications: 

•  Colin  McAuliffe,  Ryan  Karkkainen,  Chian  Yen  and  Waisman,  H.  and  Haim  Waisman  (2014), 
Numerical  modeling  of  friction  stir  welded  aluminum  joints  under  high  rate  loading.  Finite 
Elements  in  Analysis  and  Design  89:8-18. 

•  Colin  McAuliffe  and  Haim  Waisman  (2015),  A  unified  model  for  metal  failure  capturing 
shear  banding  and  fracture.  International  Journal  of  Plasticity  65:131-151. 

•  Colin  McAuliffe  and  Haim  Waisman  (2015),  On  the  importance  of  nonlinear  elastic  effects 
in  shear  band  modeling.  International  Journal  of  Plasticity  7 1 : 10-3 1 . 

•  Colin  McAuliffe  and  Haim  Waisman  (2016),  A  coupled  phase  field  shear  band  model  for 
Ductile-Brittle  transition  in  notched  plate  impact.  Computer  Methods  in  Applied  Mechanics 
and  Engineering  305,  173-195,  2016. 

•  Miguel  Arriaga  and  Haim  Waisman  (2017),  Stability  analysis  of  the  phase-field  method 
for  fracture  with  a  general  degradation  function  and  plasticity  induced  crack  generation. 
Mechanics  of  Materials,  http://dx.doi.Org/10.1016/j.mechmat.2017.04.003  . 

•  Miguel  Arriaga  and  Haim  Waisman,  Combined  stability  analysis  of  phase-field  dynamic 
fracture  and  shear  band  localization  (2017),  International  Journal  of  Plasticity  96:81-1 19. 

•  Miguel  Arriaga  and  Haim  Waisman  (2017),  A  Multidimensional  stability  analysis  of  the 
phase-field  method  for  fracture  with  a  general  degradation  function  and  energy  split.  Com¬ 
putational  Mechanics,  DOI  10.1007/s00466-017-1432-l. 

Invited  Seminars  and  Keynote  lectures  (2016-2017) 

•  Seminar:  Dynamic  Eracture  of  Metals,  College  of  Mechanics  of  Materials,  Hohai  University, 
Nanjing,  China,  July  2017. 

•  Keynote  Lecture:  Stability  analysis  of  the  Phase  Eield  method  for  fracture.  United  States  Na¬ 
tional  Congress  on  Computational  Mechanics  (USNCCM14),  Montreal,  Canada,  Jul  2017. 

•  Invited  Lecture:  A  unified  model  for  metal  failure  capturing  shear  banding  and  fracture,  lU- 
TAM  Symposium  on  Integrated  Computational  Structure-Material  Modeling  of  Deformation 
&  Eailure  Under  Extreme  Conditions,  Baltimore,  MD,  Jun  2016. 

•  Keynote  Lecture:  Stability  analysis  of  the  Phase  Eield  method  for  fracture,  European  Congress 
on  Computational  Methods  in  Applied  Sciences  and  Engineering  (ECCOMAS),  Crete  Is¬ 
land,  Greece,  Jun  2016. 
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•  Invited  Lecture:  Towards  stability  analysis  of  a  unified  model  for  metal  failure  capturing 
shear  banding  and  fracture",  lUTAM  Symposium  on  Dynamic  Instabilities  in  Solids,  Madrid, 
Spain,  May  2016. 

Conference  presentations  (2016-2017): 

•  An  overlapping  domain  decomposition  preconditioner  for  monolithic  solution  of  shearbands, 
VII  International  Conference  on  Coupled  problems  in  Science  and  Engineering,  Rhodes 
Island,  Greece,  June  2017. 

•  Stability  Analysis  of  Metals  Capturing  Brittle  and  Ductile  Fracture  through  a  Phase  Field 
Method  and  Shear  Band  Localization,  ASCE-Engineering  Mechanics  Institute  conference. 
University  of  California  San  Diego,  San  Diego,  CA,  Jun  2017. 

•  Stability  Analysis  of  Shear  bands,  ASCE-Engineering  Mechanics  Institute  conference,  Van¬ 
derbilt  University,  Nashville,  TN,  May  2016.  (Presenter:  Miguel  Arriaga) 

•  Stability  Analysis  of  Shear  bands,  Melosh  Medal  competition  in  Finite  Element  Analysis, 
Duke  University,  Durham,  NC,  April  2016.  (Presenter:  Miguel  Arriaga) 


5  Graduate  Students  Involved  in  the  project 

Two  graduate  students  were  funded  by  this  project  as  listed  below.  Both  students  have  graduated 
with  a  PhD  degree  and  one  student,  Cohn  McAuliffe,  continued  to  be  supported  by  this  project  one 
more  year  as  a  postdoc.  Furthermore,  Cohn  McAuliffe  had  a  successful  internship  with  the  Army 
Research  Lab  in  Aberdeen  Proving  Ground,  MD  during  summer  of  2012. 

Both  students  have  presented  their  work  in  domestic  and  International  conferences  such  as  the 
ASCE-Engineering  Mechanics  Institute  and  US  National  Congress  on  Computational  Mechanics. 
The  two  students  were  very  successful  in  that  their  work  was  recognized  by  the  community  winning 
several  awards  in  the  field  of  computational  mechanics.  Furthermore,  both  students  also  won 
prestigious  awards  in  the  department  for  the  best  PhD  thesis.  Please  see  below  a  detailed  award 
description. 

•  Cohn  McAuliffe  (former  Ph.D  student  and  postdoc,  now  at  Altair) 

•  Miguel  Arriaga  (former  Ph.D  student,  now  postdoc  at  Columbia  University  ME  department) 

6  Awards,  Honors  and  Appointments 

•  Haim  Waisman,  elected  chair  of  the  Computational  Mechanics  committee  in  the  Engineering 
Mechanics  Institute  of  ASCE,  September  2017. 

•  Haim  Waisman,  appointed  associate  editor  of  the  journal  of  Engineering  Mechanics,  January 
2017. 

•  Haim  Waisman,  Associate  Professor  received  tenure  at  Columbia  University,  April  2016. 
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•  Miguel  Arriaga,  won  Second  Place  in  Computational  Mechanics  Poster  Competition  at  the 
EMI-ASCE  Engineering  Mechanics  Institute  conference  at  Vanderbilt  University,  Nashville, 
TN,  May  2016. 

•  Miguel  Arriaga,  Awarded  the  Mindlin  prize  for  the  most  outstanding  PhD  thesis  in  the  Civil 
Engineering  &  Engineering  Mechanics  department  at  Columbia  University,  May  2016. 

•  Miguel  Arriaga,  Finalist  Paper  in  the  Robert  J.  Melosh  competition  for  best  paper  on  Finite 
Element  Methods,  Duke  University,  NC,  April  2016. 

•  Miguel  Arriaga,  won  Third  Place  in  Computational  Mechanics  Poster  Competition  at  the 
EMI-ASCE,  McMaster  University,  Hamilton,  Ontario,  Canada,  August  2014. 

•  Colin  McAuliffe,  Awarded  the  Mindlin  prize  for  the  most  outstanding  PhD  thesis  in  the  Civil 
Engineering  &  Engineering  Mechanics  department  in  Columbia  University,  May  2014. 

•  Haim  Waisman,  Awarded  the  EMI  Leonardo  Da  Vinci  Award  (Young  Investigator),  ASCE- 
EMI  institute,  August  2014. 
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